{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# stacking_fault_static calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The stacking_fault_static calculation style evaluates the energy of a single generalized stacking fault shift along a specified crystallographic plane.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Description updated and small changes due to iprPy version.\n",
    "- 2020-05-22: Version 0.10 update - potentials now loaded from database.\n",
    "- 2020-09-22: Calculation updated to use atomman.defect.StackingFault class. Setup and parameter definition streamlined.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- The system's dimension perpendicular to the fault plane should be large to minimize the interaction of the free surface and the stacking fault.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "First, an initial system is generated.  This is accomplished using atomman.defect.StackingFault, which\n",
    "\n",
    "1. Starts with a unit cell system.\n",
    "\n",
    "2. Generates a transformed system by rotating the unit cell such that the new\n",
    "   system's box vectors correspond to crystallographic directions, and filled\n",
    "   in with atoms to remain a perfect bulk cell when the three boundaries are\n",
    "   periodic.\n",
    "\n",
    "3. All atoms are shifted by a fractional amount of the box vectors if needed.\n",
    "\n",
    "4. A supercell system is constructed by combining multiple replicas of the\n",
    "   transformed system.\n",
    "\n",
    "5. The system is then cut by making one of the box boundaries non-periodic.  A limitation placed on the calculation is that the normal to the cut plane must correspond to one of the three Cartesian ($x$, $y$, or $z$) axes.  If true, then of the system's three box vectors ($\\vec{a}$, $\\vec{b}$, and $\\vec{c}$), two will be parallel to the plane, and the third will not.  The non-parallel box vector is called the cutboxvector, and for LAMMPS compatible systems, the following conditions can be used to check the system's compatibility:\n",
    "\n",
    "   - cutboxvector = 'c': all systems allowed.\n",
    "\n",
    "   - cutboxvector = 'b': the system's yz tilt must be zero.\n",
    "\n",
    "   - cutboxvector = 'a': the system's xy and xz tilts must be zero.\n",
    "\n",
    "A LAMMPS simulation performs an energy/force minimization on the system where the atoms are confined to only relax along the Cartesian direction normal to the cut plane.\n",
    "\n",
    "A mathematical fault plane parallel to the cut plane is defined in the middle of the system.  A generalized stacking fault system can then be created by shifting all atoms on one side of the fault plane by a vector, $\\vec{s}$.  The shifted system is then relaxed using the same confined energy/force minimization used on the non-shifted system.  The generalized stacking fault energy, $\\gamma$, can then be computed by comparing the total energy of the system, $E_{total}$, before and after $\\vec{s}$ is applied\n",
    "\n",
    "$$ \\gamma(\\vec{s}) = \\frac{E_{total}(\\vec{s}) - E_{total}(\\vec{0})}{A},$$\n",
    "\n",
    "where $A$ is the area of the fault plane, which can be computed using the two box vectors, $\\vec{a_1}$ and $\\vec{a_2}$, that are not the cutboxvector.\n",
    "\n",
    "$$A = \\left| \\vec{a_1} \\times \\vec{a_2} \\right|,$$\n",
    "\n",
    "Additionally, the relaxation normal to the glide plane is characterized using the center of mass of the atoms above and below the cut plane.  Notably, the component of the center of mass normal to the glide/cut plane is calculated for the two halves of the the system, and the difference is computed\n",
    "\n",
    "$$ \\delta = \\left<x\\right>^{+} - \\left<x\\right>^{-}.$$\n",
    "\n",
    "The relaxation normal is then taken as the change in the center of mass difference after the shift is applied.\n",
    "\n",
    "$$ \\Delta\\delta = \\delta(\\vec{s}) - \\delta(\\vec{0}).$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import shutil\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np \n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'stacking_fault_static'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "        \n",
    "# Initialize connection to library\n",
    "with open('C:/Users/lmh1/Documents/potentials_nist_gov/password.txt') as f:\n",
    "    user, pswd = f.read().strip().split()\n",
    "library = iprPy.Library(load=['lammps_potentials'], username=user, password=pswd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1 Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use.\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is loaded from the database for the prototype."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('crystal', potential=potential, family='A1--Cu--fcc', database=library)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Specify the defect parameters\n",
    "\n",
    "- __hkl__ gives the Miller (hkl) or Miller-Bravais (hkil) plane to create the free surface on.\n",
    "\n",
    "- __cutboxvector__ specifies which of the three box vectors ('a', 'b', or 'c') is to be made non-periodic to create the free surface. \n",
    "\n",
    "- __shiftindex__ can be used for complex crystals to specify different termination planes.\n",
    "\n",
    "- __a1vect_uvw, a2vect_uvw__ specify two non-parallel Miller crystal vectors within the fault plane corresponding to full planar shifts from one perfect crystal configuration to another.  \n",
    "\n",
    "- __a1, a2__ are the fractional shifts along a1vect_uvw, a2vect_uvw to apply to the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "hkl = [1, 1, 1]\n",
    "cutboxvector = 'c'\n",
    "shiftindex = 0\n",
    "a1vect_uvw = [ 0.0,-0.5, 0.5]\n",
    "a2vect_uvw = [ 0.5,-0.5, 0.0]\n",
    "a1 = 0.5\n",
    "a2 = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Modify system\n",
    "\n",
    "- __sizemults__  list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "- __minwidth__ specifies a minimum width that the system should be along the cutboxvector direction. The given sizemult in that direction will be increased if needed to ensure that the system is at least this wide.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizemults = [5, 5, 10]\n",
    "minwidth = uc.set_in_units(0.0, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.6. Specify calculation-specific run parameters\n",
    "\n",
    "- __energytolerance__ is the energy tolerance to use during the minimizations. This is unitless.\n",
    "\n",
    "- __forcetolerance__ is the force tolerance to use during the minimizations. This is in energy/length units.\n",
    "\n",
    "- __maxiterations__ is the maximum number of minimization iterations to use.\n",
    "\n",
    "- __maxevaluations__ is the maximum number of minimization evaluations to use.\n",
    "\n",
    "- __maxatommotion__ is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "energytolerance = 1e-8\n",
    "forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')\n",
    "maxiterations = 10000\n",
    "maxevaluations = 100000\n",
    "maxatommotion = uc.set_in_units(0.01, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1 sfmin.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('sfmin.template', 'w') as f:\n",
    "    f.write(\"\"\"#LAMMPS input script that performs an energy minimization\n",
    "#for a system with a stacking fault\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "<fix_cut_setforce>\n",
    "\n",
    "thermo_style custom step lx ly lz pxx pyy pzz pe\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "compute peatom all pe/atom \n",
    "\n",
    "min_modify dmax <dmax>\n",
    "\n",
    "dump dumpit all custom <maxeval> <sim_directory>*.dump id type x y z c_peatom\n",
    "dump_modify dumpit format <dump_modify_format>\n",
    "\n",
    "minimize <etol> <ftol> <maxiter> <maxeval>\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2 stackingfaultrelax()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stackingfaultrelax(lammps_command, system, potential,\n",
    "                       mpi_command=None, sim_directory=None,\n",
    "                       cutboxvector='c',\n",
    "                       etol=0.0, ftol=0.0,\n",
    "                       maxiter=10000, maxeval=100000,\n",
    "                       dmax=uc.set_in_units(0.01, 'angstrom'),\n",
    "                       lammps_date=None):\n",
    "    \"\"\"\n",
    "    Perform a stacking fault relaxation simulation for a single faultshift.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system containing a stacking fault.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    sim_directory : str, optional\n",
    "        The path to the directory to perform the simulation in.  If not\n",
    "        given, will use the current working directory.\n",
    "    cutboxvector : str, optional\n",
    "        Indicates which of the three system box vectors, 'a', 'b', or 'c', has\n",
    "        the non-periodic boundary (default is 'c').  Fault plane normal is\n",
    "        defined by the cross of the other two box vectors.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. (Default is 0.0).\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. (Default is 0.0).\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use (default is \n",
    "        10000).\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use (default is \n",
    "        100000).\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration (default is\n",
    "        0.01 Angstroms).\n",
    "    lammps_date : datetime.date or None, optional\n",
    "        The date version of the LAMMPS executable.  If None, will be identified\n",
    "        from the lammps_command (default is None).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'logfile'** (*str*) - The filename of the LAMMPS log file.\n",
    "        - **'dumpfile'** (*str*) - The filename of the LAMMPS dump file\n",
    "          of the relaxed system.\n",
    "        - **'system'** (*atomman.System*) - The relaxed system.\n",
    "        - **'E_total'** (*float*) - The total potential energy of the relaxed\n",
    "          system.\n",
    "    \n",
    "    Raises\n",
    "    ------\n",
    "    ValueError\n",
    "        For invalid cutboxvectors.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "    \n",
    "    # Give correct LAMMPS fix setforce command\n",
    "    if cutboxvector == 'a':\n",
    "        fix_cut_setforce = 'fix cut all setforce NULL 0 0'    \n",
    "    elif cutboxvector == 'b':\n",
    "        fix_cut_setforce = 'fix cut all setforce 0 NULL 0'\n",
    "    elif cutboxvector == 'c':\n",
    "        fix_cut_setforce = 'fix cut all setforce 0 0 NULL'    \n",
    "    else: \n",
    "        raise ValueError('Invalid cutboxvector')\n",
    "    \n",
    "    if sim_directory is not None:\n",
    "        # Create sim_directory if it doesn't exist\n",
    "        sim_directory = Path(sim_directory)\n",
    "        if not sim_directory.is_dir():\n",
    "            sim_directory.mkdir()\n",
    "        sim_directory = sim_directory.as_posix()+'/'\n",
    "    else:\n",
    "        # Set sim_directory if is None\n",
    "        sim_directory = ''\n",
    "    \n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "    \n",
    "    #Get lammps version date\n",
    "    if lammps_date is None:\n",
    "        lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data',\n",
    "                              f=Path(sim_directory, 'system.dat').as_posix(),\n",
    "                              potential=potential,\n",
    "                              return_pair_info=True)\n",
    "    lammps_variables['atomman_system_pair_info'] = system_info\n",
    "    lammps_variables['fix_cut_setforce'] = fix_cut_setforce\n",
    "    lammps_variables['sim_directory'] = sim_directory\n",
    "    lammps_variables['etol'] = etol\n",
    "    lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])\n",
    "    lammps_variables['maxiter'] = maxiter\n",
    "    lammps_variables['maxeval'] = maxeval\n",
    "    lammps_variables['dmax'] = uc.get_in_units(dmax, lammps_units['length'])\n",
    "    \n",
    "    # Set dump_modify format based on dump_modify_version\n",
    "    if lammps_date < datetime.date(2016, 8, 3):\n",
    "        lammps_variables['dump_modify_format'] = '\"%i %i %.13e %.13e %.13e %.13e\"'\n",
    "    else:\n",
    "        lammps_variables['dump_modify_format'] = 'float %.13e'\n",
    "    \n",
    "    # Write lammps input script\n",
    "    template_file = 'sfmin.template'\n",
    "    lammps_script = Path(sim_directory, 'sfmin.in')\n",
    "    template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables,\n",
    "                                         '<', '>'))\n",
    "    \n",
    "    # Run LAMMPS\n",
    "    output = lmp.run(lammps_command, lammps_script.as_posix(), mpi_command,\n",
    "                     logfile=Path(sim_directory, 'log.lammps').as_posix())\n",
    "    \n",
    "    # Extract output values\n",
    "    thermo = output.simulations[-1]['thermo']\n",
    "    logfile = Path(sim_directory, 'log.lammps').as_posix()\n",
    "    dumpfile = Path(sim_directory, f'{thermo.Step.values[-1]}.dump').as_posix()\n",
    "    E_total = uc.set_in_units(thermo.PotEng.values[-1],\n",
    "                              lammps_units['energy'])\n",
    "    \n",
    "    # Load relaxed system\n",
    "    system = am.load('atom_dump', dumpfile, symbols=system.symbols)\n",
    "    \n",
    "    # Return results\n",
    "    results_dict = {}\n",
    "    results_dict['logfile'] = logfile\n",
    "    results_dict['dumpfile'] = dumpfile\n",
    "    results_dict['system'] = system\n",
    "    results_dict['E_total'] = E_total\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3 stackingfault()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stackingfault(lammps_command, ucell, potential, hkl,\n",
    "                  mpi_command=None, sizemults=None, minwidth=None, even=False,\n",
    "                  a1vect_uvw=None, a2vect_uvw=None, conventional_setting='p',\n",
    "                  cutboxvector='c', faultpos_rel=None, faultpos_cart=None,\n",
    "                  a1=0.0, a2=0.0, atomshift=None, shiftindex=None,\n",
    "                  etol=0.0, ftol=0.0, maxiter=10000, maxeval=100000,\n",
    "                  dmax=uc.set_in_units(0.01, 'angstrom')):\n",
    "    \"\"\"\n",
    "    Computes the generalized stacking fault value for a single faultshift.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    ucell : atomman.System\n",
    "        The crystal unit cell to use as the basis of the stacking fault\n",
    "        configurations.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    hkl : array-like object\n",
    "        The Miller(-Bravais) crystal fault plane relative to ucell.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    sizemults : list or tuple, optional\n",
    "        The three System.supersize multipliers [a_mult, b_mult, c_mult] to use on the\n",
    "        rotated cell to build the final system. Note that the cutboxvector sizemult\n",
    "        must be an integer and not a tuple.  Default value is [1, 1, 1].\n",
    "    minwidth : float, optional\n",
    "        If given, the sizemult along the cutboxvector will be selected such that the\n",
    "        width of the resulting final system in that direction will be at least this\n",
    "        value. If both sizemults and minwidth are given, then the larger of the two\n",
    "        in the cutboxvector direction will be used. \n",
    "    even : bool, optional\n",
    "        A True value means that the sizemult for cutboxvector will be made an even\n",
    "        number by adding 1 if it is odd.  Default value is False.\n",
    "    a1vect_uvw : array-like object, optional\n",
    "        The crystal vector to use for one of the two shifting vectors.  If\n",
    "        not given, will be set to the shortest in-plane lattice vector.\n",
    "    a2vect_uvw : array-like object, optional\n",
    "        The crystal vector to use for one of the two shifting vectors.  If\n",
    "        not given, will be set to the shortest in-plane lattice vector not\n",
    "        parallel to a1vect_uvw.\n",
    "    conventional_setting : str, optional\n",
    "        Allows for rotations of a primitive unit cell to be determined from\n",
    "        (hkl) indices specified relative to a conventional unit cell.  Allowed\n",
    "        settings: 'p' for primitive (no conversion), 'f' for face-centered,\n",
    "        'i' for body-centered, and 'a', 'b', or 'c' for side-centered.  Default\n",
    "        behavior is to perform no conversion, i.e. take (hkl) relative to the\n",
    "        given ucell.\n",
    "    cutboxvector : str, optional\n",
    "        Indicates which of the three system box vectors, 'a', 'b', or 'c', to\n",
    "        cut with a non-periodic boundary (default is 'c').\n",
    "    faultpos_rel : float, optional\n",
    "        The position to place the slip plane within the system given as a\n",
    "        relative coordinate along the out-of-plane direction.  faultpos_rel\n",
    "        and faultpos_cart cannot both be given.  Default value is 0.5 if \n",
    "        faultpos_cart is also not given.\n",
    "    faultpos_cart : float, optional\n",
    "        The position to place the slip plane within the system given as a\n",
    "        Cartesian coordinate along the out-of-plane direction.  faultpos_rel\n",
    "        and faultpos_cart cannot both be given.\n",
    "    a1 : float, optional\n",
    "        The fractional coordinate to evaluate along a1vect_uvw.\n",
    "        Default value is 0.0.\n",
    "    a2 : float, optional\n",
    "        The fractional coordinate to evaluate along a2vect_uvw.\n",
    "        Default value is 0.0.\n",
    "    atomshift : array-like object, optional\n",
    "        A Cartesian vector shift to apply to all atoms.  Can be used to shift\n",
    "        atoms perpendicular to the fault plane to allow different termination\n",
    "        planes to be cut.  Cannot be given with shiftindex.\n",
    "    shiftindex : int, optional\n",
    "        Allows for selection of different termination planes based on the\n",
    "        preferred shift values determined by the underlying fault generation.\n",
    "        Cannot be given with atomshift. If neither atomshift nor shiftindex\n",
    "        given, then shiftindex will be set to 0.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. (Default is 0.0).\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. (Default is 0.0).\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use (default is \n",
    "        10000).\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use (default is \n",
    "        100000).\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration (default is\n",
    "        0.01 Angstroms).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'E_gsf'** (*float*) - The stacking fault formation energy.\n",
    "        - **'E_total_0'** (*float*) - The total potential energy of the\n",
    "          system before applying the faultshift.\n",
    "        - **'E_total_sf'** (*float*) - The total potential energy of the\n",
    "          system after applying the faultshift.\n",
    "        - **'delta_disp'** (*float*) - The change in the center of mass\n",
    "          difference between before and after applying the faultshift.\n",
    "        - **'disp_0'** (*float*) - The center of mass difference between atoms\n",
    "          above and below the fault plane in the cutboxvector direction for\n",
    "          the system before applying the faultshift.\n",
    "        - **'disp_sf'** (*float*) - The center of mass difference between \n",
    "          atoms above and below the fault plane in the cutboxvector direction\n",
    "          for the system after applying the faultshift.\n",
    "        - **'A_fault'** (*float*) - The area of the fault surface.\n",
    "        - **'dumpfile_0'** (*str*) - The name of the LAMMMPS dump file\n",
    "          associated with the relaxed system before applying the faultshift.\n",
    "        - **'dumpfile_sf'** (*str*) - The name of the LAMMMPS dump file\n",
    "          associated with the relaxed system after applying the faultshift.\n",
    "    \"\"\"\n",
    "    # Construct stacking fault configuration generator\n",
    "    gsf_gen = am.defect.StackingFault(hkl, ucell, cutboxvector=cutboxvector,\n",
    "                                      a1vect_uvw=a1vect_uvw, a2vect_uvw=a2vect_uvw,\n",
    "                                      conventional_setting=conventional_setting)\n",
    "    \n",
    "    # Check shift parameters\n",
    "    if shiftindex is not None:\n",
    "        assert atomshift is None, 'shiftindex and atomshift cannot both be given'\n",
    "        atomshift = gsf_gen.shifts[shiftindex]\n",
    "    elif atomshift is None:\n",
    "        atomshift = gsf_gen.shifts[0]\n",
    "\n",
    "    # Generate the free surface (zero-shift) configuration\n",
    "    sfsystem = gsf_gen.surface(shift=atomshift, minwidth=minwidth,\n",
    "                               sizemults=sizemults, even=even,\n",
    "                               faultpos_rel=faultpos_rel)\n",
    "\n",
    "    abovefault = gsf_gen.abovefault\n",
    "    cutindex = gsf_gen.cutindex\n",
    "    A_fault = gsf_gen.surfacearea\n",
    "\n",
    "    # Identify lammps_date version\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "\n",
    "    # Evaluate the zero shift configuration\n",
    "    zeroshift = stackingfaultrelax(lammps_command, sfsystem, potential,\n",
    "                                   mpi_command=mpi_command,\n",
    "                                   cutboxvector=cutboxvector,\n",
    "                                   etol=etol, ftol=ftol, maxiter=maxiter,\n",
    "                                   maxeval=maxeval, dmax=dmax,\n",
    "                                   lammps_date=lammps_date)\n",
    "    \n",
    "    # Extract terms\n",
    "    E_total_0 = zeroshift['E_total']\n",
    "    pos_0 = zeroshift['system'].atoms.pos\n",
    "    shutil.move('log.lammps', 'zeroshift-log.lammps')\n",
    "    shutil.move(zeroshift['dumpfile'], 'zeroshift.dump')\n",
    "\n",
    "    # Evaluate the system after shifting along the fault plane\n",
    "    sfsystem = gsf_gen.fault(a1=a1, a2=a2)\n",
    "    shifted = stackingfaultrelax(lammps_command, sfsystem, potential,\n",
    "                                 mpi_command=mpi_command,\n",
    "                                 cutboxvector=cutboxvector,\n",
    "                                 etol=etol, ftol=ftol, maxiter=maxiter,\n",
    "                                 maxeval=maxeval, dmax=dmax,\n",
    "                                 lammps_date=lammps_date)\n",
    "    \n",
    "    # Extract terms\n",
    "    E_total_sf = shifted['E_total']\n",
    "    pos_sf = shifted['system'].atoms.pos\n",
    "    shutil.move('log.lammps', 'shifted-log.lammps')\n",
    "    shutil.move(shifted['dumpfile'], 'shifted.dump')\n",
    "\n",
    "    # Compute the stacking fault energy\n",
    "    E_gsf = (E_total_sf - E_total_0) / A_fault\n",
    "    \n",
    "    # Compute the change in displacement normal to fault plane\n",
    "    disp_0 = (pos_0[abovefault, cutindex].mean()\n",
    "            - pos_0[~abovefault, cutindex].mean())\n",
    "    disp_sf = (pos_sf[abovefault, cutindex].mean()\n",
    "             - pos_sf[~abovefault, cutindex].mean())\n",
    "    delta_disp = disp_sf - disp_0\n",
    "    \n",
    "    # Return processed results\n",
    "    results = {}\n",
    "    results['E_gsf'] = E_gsf\n",
    "    results['E_total_0'] = E_total_0\n",
    "    results['E_total_sf'] = E_total_sf\n",
    "    results['delta_disp'] = delta_disp\n",
    "    results['disp_0'] = disp_0\n",
    "    results['disp_sf'] = disp_sf\n",
    "    results['A_fault'] = A_fault\n",
    "    results['dumpfile_0'] = 'zeroshift.dump'\n",
    "    results['dumpfile_sf'] = 'shifted.dump'\n",
    "    \n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = stackingfault(lammps_command, ucell, potential, hkl,\n",
    "                             mpi_command = mpi_command,\n",
    "                             sizemults = sizemults,\n",
    "                             minwidth = minwidth,\n",
    "                             a1vect_uvw = a1vect_uvw,\n",
    "                             a2vect_uvw = a2vect_uvw,\n",
    "                             cutboxvector = cutboxvector,\n",
    "                             shiftindex = shiftindex,\n",
    "                             a1 = a1,\n",
    "                             a2 = a2,\n",
    "                             etol = energytolerance,\n",
    "                             ftol = forcetolerance,\n",
    "                             maxiter = maxiterations,\n",
    "                             maxeval = maxevaluations,\n",
    "                             dmax = maxatommotion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['E_gsf', 'E_total_0', 'E_total_sf', 'delta_disp', 'disp_0', 'disp_sf', 'A_fault', 'dumpfile_0', 'dumpfile_sf'])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1 Define units for outputting values\n",
    "\n",
    "- __length_unit__ is the unit of area to display delta displacemets in.\n",
    "\n",
    "- __area_unit__ is the unit of area to display fault area in.\n",
    "\n",
    "- __energyperarea_unit__ is the energy per area to report the surface energy in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "length_unit = 'nm'\n",
    "area_unit = 'nm^2'\n",
    "energyperarea_unit = 'mJ/m^2'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2 Print $A_{fault}$, $E_{gsf}$, and $\\Delta\\delta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Values for fractional shift = (0.500000, 0.000000)\n",
      "A_fault =     5.365198866947918 nm^2\n",
      "E_gsf =       924.3579463239225 mJ/m^2\n",
      "delta_disp =  0.015655877018008368 nm\n"
     ]
    }
   ],
   "source": [
    "print('Values for fractional shift = (%f, %f)' % (a1, a2))\n",
    "print('A_fault =    ', uc.get_in_units(results_dict['A_fault'], area_unit), area_unit)\n",
    "print('E_gsf =      ', uc.get_in_units(results_dict['E_gsf'], energyperarea_unit), energyperarea_unit)\n",
    "print('delta_disp = ', uc.get_in_units(results_dict['delta_disp'], length_unit), length_unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
